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Dynamical evolution of rotating dense stellar systems with 
embedded black holes 
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ABSTRACT 

Evolution of sclf-gravitating rotating dense stellar systems (e.g. globular clusters, 
galactic nuclei) with embedded black holes is investigated. The interaction between 
the black hole and stellar component in differential rotating flattened systems is fol- 
lowed. The interplay between velocity diffusion due to relaxation and black hole star 
accretion is investigated together with cluster rotation using 2D+f Fokkcr-Planck 
numerical methods. The models can reproduce the Bahcall-Wolf solution / oc E 1 / 4 
(n oc r~ 7 / 4 ) inside the zone of influence of the black hole. Gravo-gyro and gravother- 
mal instabilities conduce the system to a faster evolution leading to shorter collapse 
times with respect to the non-rotating systems. Angular momentum transport and 
star accretion support the development of central rotation in relaxation time scales. 
We explore system dissolution due to mass-loss in the presence of an external tidal 
field (e.g. globular clusters in galaxies). 
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1 INTRODUCTION 

Observational analysis of globular clusters (GCs) has been 
considerably improved in the last year s thanks the Hub- 
ble Space Telescope (HST), (cf. e . g. IPiotto et al.l l2002i ; 
iRich et alj|2005l ; iBeccari et afl|2006l : iGeorgiev et"ajTl2009ft . 
They have been used to obtain luminosity functions and 
derived mass functions, color-magnitude diagrams (CMDs) 
and population and kinematical analysis leading to a bet- 
ter understanding of their evolutionary processes. At the 
same time, new questions have been opened due to the, 
in the past unimpressive, complexity of this stellar sys- 
tems. Nevertheless, analysis of rotation in these systems 
has been rarely carried on. On the other side, there are 
observational evidences for the existence of intermediate- 
mass black holes (IMBHs) in GCs, although their origin 
is as yet not very clear. Local core collapsed GCs are ex- 
pected to harbor IMBHs due to their high central densi- 
ties, but regarding its detection, it was argued that 'non- 
collapsed' projected density profiles which evolve harbor- 
ing IMBHs fit well by medium-concentra t ion King m od- 
els (|Baumgardt et al.ll2005r i. iGerssen et all (|2002l , |2003h re- 
ported the kinematical study (based on HST spectra) of 
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the central part of the collapsed GC M15. They proposed 
the presence of an IMBH (M bh = 3.9 ■ 10 3 M Q ) in the cen- 
tral region of M15. A single or binary IMBH could ac- 
count for the net rotation observed in the center of M15 
dGebhaxdt et al.ll2000l:lGerssen et aI|2002i;lMuler fc Colbert! 
12004 iKiselev et all |200S|) . iMaccarone fc ServillatJ (|2008h 
shows the possible presence of an IMBH in NGC 2 808 with 
a mass of ~ 2.7 x 10 3 M Q and iNovola et all (120081 ) have re- 
ported ab out the BH in omegaCen with a mass of ~ 1O 4 M . 
However. iBaumgardt et all (|2003T ) have shown through self- 
consistent A-body computations treating stellar evolution 
and with a realistic IMF, that the core collapse profile of a 
star cluster with an unseen concentration of neutron stars 
and heavy-mass white dwarfs can explain the observed cen- 
tral r ise of the mass-to-light ratio (see also iMcNamara et ail 
l2003h . Similarly, a dense concentration of compact rem- 
nants might also be responsible for the high mass-to-light 
ratio of the central region of NGC 6752 s een in pulsar tim- 
ings (|Ferraro et"al]|2003! ; [Colpi et al.ll2003l). Out s ide ou r own 
galaxy, iGebhardt et all (|2002l , 120051 ): IZahariias! (|2008l ) have 
reported evidence for a 20000 M BH in the M31 globular 
cluster Gl. 

Chandra and XMM-Newton observations of ultralumi- 
nous X-ray sources (ULXR) give also evidence for the ex- 
istence of IMBHs in dense star clusters, which are often 
associated with young star clusters and whose high X-ray 
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luminosities in m any cases suggest a compact object mas s 
of at least 10 2 M o (jEbisuzaki etal.ll200ll ; [Miller et al.ll2003l ). 
Furthermore, some of the ULXR sources detected in other 
galax ies are may be accreting IMBHs fe.g.. lMiller fe Colbert! 
120041 '). altho ugh the majority could be likely stella r-mass 
black holes (|King et al.ll200ll ; iRappaport et al.ll2005l ). 

On the other side, the centers of most galaxies embed 
massive BHs. It has been evidenced by HST measurements 
in the last years, and since theoretical modeling of measured 
motions request for the presence of cen tral compact dark ob- 
jects with a mas s of ~ 10 6 to 10 9 M B (iFerrarese et al Il200ll ; 
iGebhardtl |2002| ; iPinknev et all 120031 : kormendv, J. 1 12004 ). 
Ground-based IR observations of the fast orbital motions 
of a few stars in the Milky Way have l ead to the detection 
of a 3 - 4 x 10 6 M B BH in its cent er (|Schodel et al.ll2003l ; 
iGhez et al.ll2004l ; lEckart et al.ll2004T ). Moreover, BH demo- 
graphics have lead to correlations between the BH mass 
and the luminosity of its ho st bulge or elliptical galaxy 
|Kormendv fc Richstonelll995l ). and between BH mass and 
the velocity dispersion of its host bulge, as Mbh oc a a 
(|Ferrarese fc Merritt! 120001 ) . This leads to a strong link be- 
tween BH formation and the properties of the stellar bulge, 
like the formation of density cus ps (Bahcall-Wolf solution), 
which have been investigated bv lSchodel. Merritt fc Eckartl 
(200I) 

Dynamical modeling of collisional stellar systems (like 
galactic nuclei, rich open clusters, and rich galaxy clusters) 
still possess a considerable challenge for both theory and 
computational requirements (in hardware and software). On 
the theoretical side the validity of certain assumptions used 
in statistical modeling based on the Fokker-Planck (FP) 
and other approximations has not been fully investigated. 
Stochastic noise in a discrete A-body system and the impos- 
sibility to directly model realistic particle numbers with the 
presently available hardware, are a considerable challe nge 
for the computational side (but see lBerczik et all (|2006t )) 

While all work known to the authors at this moment 
concentrates on self-gravitating star clusters, the improve- 
ment of our knowledge and methods in the field of rotat- 
ing dense stellar systems is extremely important for galac- 
tic nuclei, too, where a central star-accreting black hole 
comes i nto the game (some stationary modeling exists , 
such as iDuncan fc Shapiro! I1983L lOuinlan fc Shapiro! 1 1990l . 



iMurphv B. et all Il99ll . IFreitag fc Bend I2002T ). Direct inte- 
gration of or bits (A-Body method) has been applied t o 
the problem (|Giiltekin et alj|2004l ; iBaumgardt et al.ll2005t) . 
However, A-Body simulations only provide a very limited 
number of case studies, due to the enormous computing 
time needed even on the GRAPE computers. Moreover, re- 
cent investigations show that in young dense clusters, super- 
massive stars may form through runaway merging of main- 
sequence stars via direct physical collisions, which may then 
collapse to form an IMBH. The collision rate will be greatly 
enhanced if massive stars hav e time to reach the core be- 
fore exploding as supernovae |portegies Zwart et al.l |2004| ; 
iGtirkan et "aI1l2006l ; IFreitag et al.ll2006l ). Therefore it is very 
urgent to develop reliable approximate models of rotating 
star clusters with black hole, which is subject of the present 
work. 

A 2D FP model has been wor ked out for the cas e of ax- 
isymmetric rotating star clusters (|Paper itlPaper 11) . Here, 
the distribution function is assumed to be a function of en- 



ergy E and the ^-component of angular momentum (J z ) 
only; a possible dependence of the distribution function on 
a third integral is neglected. As in the spherically symmetric 
case the neglect of an integral of motion is equivalent to the 
assumption of isotropy, here between the velocity dispersions 
in the meridional plane (zu and z directions); anisotropy be- 
tween velocity dispersion in the meridional plane and that 
in the equatorial plane ((^-direction), however, is included. 

We realize that the evolutionary models provided by 
us for rotating dense stellar systems are difficult to use for 
direct comparisons with observations, because they are not 
easily analytically describable. But they are the only ones 
which fully cope with all observational data available nowa- 
days (full 3D velocity data, including velocity dispersions in 
zu and (^-direction, rotational velocity, density, all as full 2D 
functions of zu and z, see Fiestas et al. 2006). No other evo- 
lutionary model exists so far which is able to provide this 
information. With t he advent of our new post-collapse and 
multi-mass models |Paper III; iPaper IlJ) and the inclusion 
of stellar evolution and binaries (work in progress) we will 
be able to deliver even more interesting results. Already the 
existing A-body study (|Ardi. Spurzem fc Mineshige] [2005) 
shows that rotation not only accelerates the collisional evo- 
lution (but see Ernst et al. 2007) but also leads to an in- 
creasing binary activity in the system. 

A description of the method used in the present work 
is made in Section [21 Section describes the initial con- 
figurations of the models and numerical tests of the code, 
Section [4] presents the main results in the isolated case first 
reproducing the spherically symmetric model (no rotation), 
and secondly, giving a description of rotational behavior of 
axisymmetric systems in relaxation time scales, emphasizing 
the interplay between the dynamical evolutionary processes. 
Section [5] explores the system dissolution due to the tidal 
field of a parent galaxy. Section [6] gives the conclusions and 
further plans. 



2 THEORETICAL MODEL 
2.1 Equations and assumptions 

The pioneering work of lGoodmanl (| 19831 ) , in his unpublished 
thesis, and the fur ther development o f the Fokker-P lanck 
method made by 



lEinsel fc Spurzem! dl999. Paper J ) and 
iKim et al.l (|2002 Paper III 12004 Paper lit , have brought 
the treatment of the axisymmetric rotating case to a 
newly state of interest, which follows the evolution of self- 
gravitating rotating systems driven by relaxation effects and 
its consequences for the stellar redistribution and shape of 
the system. 

Evolution of the distribution function /(r, v) Q of stars 
in phase space (r, v) under the influence of the potential (j)(f) 
is described by the Boltzmann-equation 



^+^V r +--V v / = (^-)coll 

at m at 



(1) 



with space and velocity coordinates, r and v respectively. 
The force F = —mV T (f> is applied on stars of mass m. The 



a proper definition of / corresponding to the initial conditions 
for this study is given in Eq. 1291 
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term on the right side of Eq.[T]takes into account the changes 
in / due to collisions (not real collisions but stellar scatter- 
ings which cause deviations in the orbits) . The collision term 
is given through the -local- Fokker-Planck approximation 



21) = _^_ (/(Au M >) + i^L_ 

dt /con dv» UK " 2dv»dv 



(/(At^Au")) (2) 



where fi = 1,2,3 and v = 1,2,3 (tensor notation). w M gives the 
velocity in Cartesian coordinates. The first order diffusion 
coefficients (Aw M ) describe the dynamical friction, while the 
second order ones (Aii^Au") give the real velocity diffusion. 

In order to obtain the solution of the Fokker-Planck 
equation following assumptions are taken into account: 

• cluster evolution time scales are in following relation: 

tdyn <S trh <S tcl (3) 

where t c i represents the cluster age. t r h is the half-mass re- 
laxation time, following ISpitzer fc Hart! |l97ll ): 



trh 



0.138, 



f Nrl 



(4) 



N is the number of particles (stars), G the gravitational 
constant, m the mean stellar mass, InA is the Coulomb log- 
arithm (A = 0.47V is used here) and r n the half-mass radius. 
The dynamical time is given by 



tdyn — 1-58 



(5) 



GM 

where M is the total mass of the cluster. 
The system evolves slowly through diffusion in a sequence of 
virtual equilibrium states. In a time t r h information about 
the initial configuration is lost due to relaxation. A prove of 
this statement is shown in Section f3. 21 

• the solution is given for small-angle scatterings 
(Av/v <C 1), i.e., for changes of v to v + Av. 

• there is no correlation between collisions (like in three- 
body collisions), which could be important for the energy 
generation in the core that can reverse the collapse. 

• no binaries and stellar evolution are considered. Thus, 
binary heating due to 3-body encounters is neglected . 
Note that binary heating can reverse coll apse (|Hutlll985l : 
iMcMillan. Hut fc Makinoill990l ; iKim et al.ll2002 Paper it) 

• We neglect any recoil of the BH as a result of accretion 
and three-body interactions. If these were to be taken into 
account the initial BH seed would end up kicked out of the 
system well before it has a chance to significantly grow in 
mass. Realistic effects on the dynamical evolution are being 
studied by using NBody realisations of the systems, to be 
published in a forthcoming paper. 

• the initial BH mass (MbhJ is much smaller than the 
cluster mass M c i 

• The distribution of stars is represented by an equal- 
mass particle system, which is initially axisymmetric in 
space and is able to develop anisotropy in velocity space. 
No stellar spectrum is included in this model, with the aim 
to test the model without large complexity. 

The classical isolating integrals of a general axisymmet- 
ric potential <j>, in cylindrical coordinates (zu, z), are the en- 
ergy per unit mass: 

1 



where <f) c \(zu,z) is the potential of the stellar system and 
</>bh(s7, z) = —GAIbh/r is the BH-potential (r 2 = zu 2 + z 2 )\ 
and the component of angular momentum along the z-axis 
per unit mass, given by 

J z — zuve v (7) 

v v = ve v is the velocity component in azimuthal direction. 
E and 4> are negative for all particles. 

Conservation of E and J z is used in the solution of the 
Fokker-Planck equation, which becomes a non-linear second 
order integro-differential equation (the diffusion coefficients 
of Eq.[5]are expressed in terms of integrals over the local field 
star velocity distribution function). This integrals are given 
by the Rosenbluth potentials (Rosenbluth et al., 1987). A 
derivation of th e diffusion coefficients in terms of E and J z 
can be found in lEinsel fc SpurzemI (|l999. Paper Jl . 

In axisymmetric systems, although / can be approxi- 
mately representable as a function of E, J z and t (except 
for very special forms of the potential), numerical evidence 
demonstrates that axisymmetric potentials can support or- 
bits which have three integrals of motion: E, J z and a 
third integral commonly designated J3. That is, the typi- 
cal orbit does not spread uniformly over the hypersurfacc 
in phase space defined by its E and J z but is confined to 
lower-dimensional subset ('non-ergodic' orbits on their EJ Z 
surfaces). A solution of the orbit-averaged FP equation in 
energy-momentum space may represent an artificial case of 
a true point-mass system, since in the axisymmetric poten- 
tial a third int egral of motion c ould restrict particle motion 
in phase space (|Goodmanlli983T ). Since on one side, the inner 
parts of the cluster are dominated by relaxation effects and 
the third integral can be neglected, due to the efficiency of 
diffusion in these regions; on the other side, the outer region 
of the cluster can be strong influenced by the third integral, 
as radially biased anisotropy dominates this region. 

In the present study, non-ergodicity on the hypersur- 
face (given by E and J z ) is neglected due to any third 
integral I3. The potential close to the BH is spherically 
symmetric (~ 1/r), and 73 could be fairly approximated 
by J 2 , since less radial orbits are expected in this regions 
llAmaro-Seoane. Freitag fc Spurzem 2004: B aumgardt et all 
l2004h . which are preferentially disrupted by the BH. The 
angular momentum J z is here good represented by its max- 
imum value (J™ ax ). Nevertheless, possible existing merid- 
ional circular orbits can not be distinguished by our model 
and would be treated as radial orbits (for example for their 
accretion). 

In terms of integrals of motion, the Boltzmann equation 
(Eq. [1} is expressed in the axisymmetric system as: 

21 + 2±21 = (21\ 
dt dt dE \at /con 



(8) 



the dependence on J z is given implicitly by (f>\ and the colli- 
sional term of Eq. [2] can be expressed in terms of E and J z 



(21) = i[. 

Vat /coll VI 



-^({AE)Vf)-^-((AJ.)Vf) + 



1 d 2 



+ ^^L((( AE f)Vf) + ^({AEAJz)Vf)+ 



E = ^v 2 + </> cl (zu, z) + 0bh(^, 2) 



(6) 



2dE 2 
1 d 2 



2<9J Z 



({(AJ z Y)Vf) 



(9) 
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with the volume element in velocity space given by V = —. 

The vast dynamical parameter range o f relaxed and 
unrela xed cluster systems (as presented in iFiestas et al.l 
(120061 )1 was specially treated applying appropriate bound- 
ary conditions at the inner potential cusp of the BH and 
the outer cluster tidal boundary (in the presence of a parent 
galaxy). A double-logarithmic (zu, z) space grid is used, and 
the Fokker-Planck equation is written in a dimensionless flux 
form by introducing the dimensionless energy 



X{E) = ln( 



E 



24> c -E Q -E 



(10) 



Eo is a characteristic energy, which allows a higher resolution 
at higher energy (and low angular momentum) levels, as 
well as in the outer parts of the system (halo), where the 
proportionality X — ln\E\ improves the spacing of the radii of 
circular orbits with given energies in the direction of the tidal 
boundary. The dimensionless angular momentum is given by 



Y{J„E) 



Jz 



u z 



(11) 



Each time step r c i rc (E) and J™ ax (i?) are determined in 
the equatorial plane from the evolving potential by a simple 
Newton-Raphson scheme, using the relations: 



(E - <f}(zu cilc , 2 = 0)) = T-OTcircT^- 

2 ozu 



(12) 



in order to get rzr c i rc (or r cr ic, due to z = 0), and computing 
Jr*(E), using: 

(J z max (£)) 2 =r c 3 irc # (13) 



2.2 Diffusion and loss-cone accretion 

Given the values of E and J z , the orbit average of the Fokker- 
Planck equation in the form of Eq.|8]is obtained by integrat- 
ing it over an area P(E, J z ,t) of the hypersurface in phase 
space, given by 



P(E,J z ,t) = 4tt^ 



dujdz 



(14) 



This weighting factor gives also the number of stars in the 
system taking part in the diffusion, as 



N(E,J.,t) = P(E,J.,t)-f(E,J z ,t) 



(15) 



A(E, J z ) is given by the intersection of the hypersurface with 
the roz-plane, where the sum of the squares of the velocity 
components are non-negative: 

A(E,J Z )= Uwz) | \vl+vl = E-^-^^o) (16) 

The condition (|16p is rastered numerically in the code by 
given E and J z . 

In a general axisymmetric potential, almost none of the 
orbits are closed, so that the orbital period is not well de- 
fined. There exist two different epicycle periods, one each 
for oscillations in the w and z directions. The orbit aver- 
age is taken over a time that is larger than both and is the 
time required for the orbit to spread uniformly over the area 
A(E, J z ), only because of encounters, i.e. in a relaxation time 
scale (if the third integral is well conserved). 



The Fokker-Planck equation is solved numerically in 
flux conservation form, following : 

dFy 
dY ' 



d£ = 1 dFx 
dt p [ dX 



(17) 



p is the phase volume per unit X and Y , with particle flux 
components in the X and Y directions: 



Fx = -D xx -^-D XY ^-D x f; 
Fy = -D YY ^-D Yx ^-D Y f 



(18) 
(19) 



The orbit-averaged flux coefficients Da are derived from the 
local diffusion coefficients and transfor med to dimension- 
less variables -Dx, Dy, Oxx, Dyy, Dxy (|Einsel fc Spurzeml 
1 1999. Paperl ). 

The loss-cone limit is defined by the minimum angular 
momentum for an orbit of energy E: 



JT m (E) = r d ^2(E - GM bh /r d ) (20) 

where r d is the disruption r adius of the BH, calculated fol- 
lowing iFrank fe Reesf |l976h : 

r d ocr,(A/ bh /m,) 1/3 (21) 

r* and m* are the stellar radius and mass, respectively. Their 
adopted values are given in Section \3. II 

The central potential cusp of an embedded massive BH 
disturbs the redistribution of stars due to collisional inter- 
actions. Thus, following assumptions are taken in order to 
develop the structural parameters of the cluster: 

• a seed initial BH mass, which is much larger than a stel- 
lar mass, is calculated numerically using a first perturbation 
of the potential in the initial models (see also Section 13.11) 

• accretion is driven by angular momentum diffusion. A 
star is completely accreted if its z-component of angular 
momentum is less than J™ m , which defines the loss-cone 
boundary. Energy diffusion for accretion is neglected be- 
cause the changes in E due to collisions are considered 
small in comparison to the changes in angular momentum 
(|Cohn fc Kulsrudlll978l ). 

• the distribution function vanishes for J z > J™ ax and 
J z < -J z max . 

• the central BH grows slowly through accretion of stars 
leading to a new distribution f(E, J z ) and a new cj>(w, z). 

• The unit of time is proportional to the relaxation time 
at the influence radius of the BH r a , defined as the radius, 
where the mass of the cluster equals Mbh- The time step is 
given by At = f [t)r T& , where: 



0.338cr 3 



n a (Gm*) 2 ZnA 



(22) 



<7 a and n a are the velocity dispersion and density evaluated 
at r a . ^(O) depends on the initial model and is increased from 
time to time by a factor of 4/3 in order to have a fractional 
increasing of central density of between 2 and 4 % per time 
step. In analogy with the computations of lCohnl (|l979l ). one 
Vlasov step, i.e. one recomputation of the potential, follows 
every Fokker-Planck (diffusion) step. 

A schematic diagram of the numerical XV-grid, as used 
for the solution method of the discretized Fokker-Planck 
equation, is shown in Fig. [1] Energy limits are the central 
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Figure 1. Schematic diagram of the numerical XY-grid and def- 
inition of the loss-cone. Half of the grid is shown in the left side 
of the figure. The other lower half corresponds to negative val- 
ues of Y and is symmetric with respect to the Y = 0-axis. The 
dark shaded area represents the loss-cone in X — Y-space, and is 
limited by Ymi n . Stars are able to go into and out of it through 
angular momentum diffusion, as shown in the right side of the 
figure, for one grid cell. 



potential (X(<f> c )) and the tidal energy (X(Etid)), angular 
momentum limits are the maximum values of Y in both di- 
rections (Y = ±1). For the purpose of this illustration, only 
the upper half of the grid is shown (0 ^ Y < +1), since 
the lower half (— 1 < Y < 0) is symmetric with respect to 
the Y — 0-axis. Angular momentum diffusion of stars into 
and out of neighboring cells is illustrated in the right part of 
the figure. The distribution function is defined in the center 
of each cell, while the diffusion terms are computed at the 
cell boundaries. As shown in Fig. Q]the fluxes at the Y=l 
boundaries are set to zero. Open boundaries are the loss- 
cone region and E = _B t id in the tidally limited models. The 
loss-cone is limited by Y mm (.E) and the limit of maximum 
energy (E((f> c )). The first derivative of f with respect to E 
is non-zero at the boundary of each grid cell and is evalu- 
ated just inside the boundary in order to obtain accurate 
escape fluxes. During the solution method, the whole grid is 
rastered and the angular momentum fluxes are saved each 
time step. In case of isolated models X(E t id) is set to 1/10 of 
the potential at the tidal radius of the corresponding King 
model. Evaporation of stars in isolated systems is neglected. 
In tidally limited systems, stars are able to escape the sys- 
tem through the tidal boundary at X(E t id), influenced by 
the potential field of the parent galaxy. 

The flux term Fy , per unit energy and unit time across 
the angular momentum boundary, used to compute contri- 
bution to stellar accretion, is given by the second order an- 
gular momentum diffusion term in Eq. 1191 



Fy 



-D y 



BY 



(23) 



The dimensionless change in f(X,Y), due to diffusion 
in the inner/outer direction is then given in a discretized 
form by: 

A/ _ At AF Y 1 

T ~ ~P ay/ 



(24) 



where At is the time step, P(X, Y) the phase space volume 
and AFV is the neto angular momentum flux in each cell. 
Energy fluxes are negl ected for accretion since they are small 
in comparison to Fy (|Cohn fc Kulsrudlll97ct ). 

After redistribution of orbits, due to small-angle colli- 
sions, those with Y ^ Ymin lie in the loss-cone. Using the 
time scales of replenishme nt tj n oc (y mm ) 2 ffl and loss -cone 
depletion t out oc (Y diff ) 2 (|Lightman fc Shapiro! ll977T ), the 
contribution to f(X, Y) of accretion should take into ac- 
count the ratio 



tout (y diff ) 2 

tin ^ (^ min ) 2 



(25) 



ymm _ jmin^jmax ana i wg (J eno t e ^jjg dimensionless angular 
momentum diffusion term due to gravitational scattering per 
time step.(Y diff ) 2 as < (AY) 2 > 

If q < 1, most loss-cone stars remain inside and 
Af(X, Y) is well represented by the flux of Eq. [24] but, be- 
cause stars could be scattered out of the loss-cone to orbits of 
Y > Ymin in an orbital time scale, a correction to the angular 
momentum flux is necessary. In the classical approximation, 
accretion of stars inside Y m i n leads to an 'empty loss-cone' 
{f(X, Y) = 0), but this is not a realistic boundary condition. 
If g > 1, the angular momentum diffusion term is larger than 
the loss-cone opening, so that most stars manage to scatter 
out of it and are not accreted. q(E CT it) = 1, defines a criti- 
cal energy (i?crit) at a radius r cr i t , in the equatorial plane, 
which marks the transition between the 'full' and 'empty' 
loss-cone regimes. 

This correction is implemented in the code by using a 
probability of accretion P a (X, Y, Y dlff ), that a star do not 
escape from the loss-cone due to angular momentum change 
once it is inside, or it can fall into the loss-cone from outside. 
In the latter probability of accretion 1 — P a is applied. 

The probability that an orbit with dimensionless angular 
momentum Y suffers a change AY < \Y — Y m i n |, where 
| Y — Ymin| i s the distance to the loss-cone boundary, is given 
by: 



Pa(Y) 



Y-Y min \ 



7^d«F exp( - 



Y' ; 



(ydiff): 



;)dY' 



(26) 



i.e., centered at each (X, Y) grid cell a Gaussian distribution 
of orbits in Y with the dispersion Y dlff is assumed. Finally, 
the contribution of f(X, Y) to accretion, at each energy- 
angular momentum grid cell is given by: 



Pa(Y)(f° U + A/) 



(27) 



yoid _|_ A j j g £ ne resu it f redistribution of stars due to 
diffusion processes. In order to get the total BH accretion 
mass, the distribution function f(E, J z ) is integrated in real 
and phase space, as: 



^tid „tid 

zudvj j dz 



7* 

Jo 



t2it r tid r +J * i 

I — / dE dJ* A/ acc (E,J z) J 



(28) 



2 the square root dependence of y mm on time reflects the fact, 
that entry into the loss-cone is a diffusive process 
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Figure 2. /( J z ), at constant energy E, for an initial model (Wo = 
0.6, loo = 0.9). From top to bottom the curves have smaller \E\, 
listed on the left column. Grid dimension is (200 X 201) in (E, J z )- 
space used to construct the models. 



where ro tld and z tld give the tidal cluster radius in zxj and 
z directions respectively. A factor of 2 x 2tt = 4tt is due to 
consideration of positive and negative zenithal coordinates 
and that the azimuthal component is symmetric. More- 
over, J™ ax = zuy/2(E — </>). The accretion mass is added to 
M bh = Mbh d + Am acc and furthermore M cl = M°, ld - Am acc . 



3 NUMERICAL RESULTS 
3.1 Initial conditions 

As initial configurations, truncated King models with added 
bulk motion are used. Their adopted distribution function 
is 



f(E,J z 



exp 



, -Etid — E . 
exp( = ) 



E < E m 



f(E,J z )=0 



, E > Etid 
(29) 

where cr c is the central ID velocity dispersion and fio is 
an angular velocity. In Fig. [2] f(J z ) at constant energy E 
against J z , for an initial high rotating model (Wo = 0.6, 
ujq = 0.9), is plotted. f(J z ) covers a wide range of values 
in logarithmic scale, and J z varies from negative to positive 
values, accordingly to two directions of rotation around the 
z-axis. J z = represents stars on radial orbits, while J™ ax , 
circular orbits. Note that the angular velocity flo in Eq. [29] 
is given by the slope of / in each curve of constant energy 
(a property of King models). The isoenergy sections become 
shorter, due to the smaller possible J™ 3 * at higher absolute 
values of energy. 

The system of units for the initial King models is given 

by 



M cli 



(30) 



where M c n is the initial mass of the cluster and r c i the initial 
core (King) radius 



9a c 



(31) 



Where n c is the central number density. 

The initial conditions of each model are given by the 
pair (Wo,uo), see Table [T] Here, Wo is the familiar King 
parameter 



Wo 



0(7tid) - 0c) 



(32) 



where 4 > { r tid) is the potential at the cluster tidal boundary 
and <fic the central potential. 



LOO = 



(33) 



is the initial rotational parameter. Radii are given in units 
of the initial cluster core radius. Table [1] shows initial pa- 
rameters of the models. Intermediate models (Wo = 6.0) are 
expected to reproduce current evolutionary states of most 
GCs. Their initial concentrations decrease and their dy- 
namical ellipticities increase, the higher t he initial ro t ation . 

b/a [f|are calculated following Goodman! |l983) 



('A 



= 1 



as defined in lEinsel fc Spurzeml (|l999. Paper I ). Tidally lim- 
ited models presented in Section [5] are denoted by MIT to 
M5T. 

There are two scaling parameters in the simulation: (i) 
the particle number N, which defines the mass of the single 
star to the total mass of the system, i.e. m = M c i.JN\ (with 
M c ij = 1 in our units) and (ii) the mass of the black hole 
with respect to the total mass, which we specify in terms 
of initial seed black hole mass j3 = M-^JM C \ = 1 x 10 -5 . 
Since we aim to study the dependence of the standard loss 
cone accretion model on the new physical situation of a ro- 
tating axisymmetric system surrounding the black hole, we 
fix /3 and make our models scale-invariant to the particle 
number as long as the pure point-mass interactions are con- 
sidered. The choice of j3 is somewhat arbitrary, and it does 
not c hange much the physics, because Fokker-P lanck models 
(e.g. lAmaro-Seoane. Freitag fc Spurzeml (|2004l )) show that 
the time evolution of the black hole and star cluster do not 
depend sensitively on the initial seed. In order to prove this 
statement, we performed tests with different values of /3, as 
shown in Fig. [3] All tests show a common final Mbh and the 
disruption rates follow the same self-similar evolution. For 
comparison, we include physical units in the right Y-axis 
and top X-axis. Here we applied our models to a massive 
globular cluster by using Mi = 5 x 10 6 M Q . The initial t rb 
was calculated following Eq. 2J after scaling the half-mass 
radius to the initial King radius using Table [T](rh = 2.7pc). 
In this equation, we use iVj = 5 x 10 6 and thus, m = 1M©. 
We obtain t r h ~ 1.4 x W 9 yr, which is a typical value of 
massive galactic globular clusters. In our models the ini- 
tial seed represents technically a small perturbation in the 
central potential of our initial model, which accretes mass 
corresponding to the loss-cone accretion onto a fixed black 
hole, scaled down to our system (see below). This unphysical 
assumption for the initial accretion is used, because our goal 
is to study the long-term evolution of the system (which ap- 
proaches a self-similar solution) and not the initial growth 
process. 

For stellar disruption we do have another parameter 
though, which is the stellar radius in our simulation units. 
This radius, defines then the disruption radius through 



axis ratio of an oblate spheroid 
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Table 1. Parameters of initial models used in the simulations. 
Column 1: model identification name; Column 2: King potential; 
Column 3: dimcnsionlcss rotation; Column 4: concentration; Col- 
umn 5: ln(r I1 /r c i); Column 6: dynamical ellipticity; Column 7: 
initial half-mass relaxation time in code units; Column 8: rate of 
rotational energy to kinetic energy. 



Figure 3. Evolution of model Ml (6.0, 0.0) for different values of 
/3 (a = 2 X 10 — 8 ). (a) BH mass against time. Left Y-axis shows 
code units (M/M c n) , right Y-axis transforms them to units of 
Mq for a M cl . = 5 X 10 6 Mq. (b) Evolution of disruption rates. 
Left Y-axis shows code units (d(Mt,h/-Mcii)/ii(t/*rhi))- Right Y- 
axis transforms them to units of Mq/j/t\ 
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Figure 4. Evolution of model Ml (6.0, 0.0) for different values of 
a (/? = 1 X 10 -5 ). (a) BH mass against time, (b) Evolution of 
disruption rate. Axis are labelled as in Fig. [3] 



Eq. [21] which grows in time due to black hole growth. We 
use as here for most simulations a value a = r*/r c i = 
2 ■ 1(T 8 . In the following we use fiducial values of (a, f3) 
= (2 X 1CT 8 , 1 x 10~ 5 ) for most of our models, which (taking 
r* as solar radius) define parameters for a globular cluster. 
We have tested the variation of a, the results are shown in 
Fig. [4] Changes of a by one order of magnitude influence 
our observables Mbh and dM^/dt only by a factor of 2-3. 
That suggests that our result can be applied by scaling to a 
wider range of astrophysical systems including galactic nu- 
clei. Moreover, for a direct application to globular clusters 
our initial model is unphysical, because the seed black hole 
is not fixed, and it may grow or be ejected by close three- 
body encounters, all effects we are currently not taking into 
account. However, no matter what is the growth mechanism, 
if the black hole remains in the cluster and grows, we think 
that our scale-invariant solution should be reached. 
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Table 2. Convergence analysis of relative errors performed by 
the code. Column 1: Grid size; Column 2: Relative energy er- 
ror; Column 3: Relative mass error; Column 4: Relative angular 
momentum error 

3.2 Numerical tests 

We use a grid size of N x = 200, N Y = 201, iV w = N z = 200 
obtaining errors of angular momentum, mass and energy as 
shown in Table [5] for a typical model (our reference model, 
Wo — 6) by the time the central density had decreased by 
about 2-3 orders of magnitude, during core expansion. As 
can be read from Table [5] a grid convergence study shows 
that shorter grids are not sufficient to keep the errors small 
and in order to achieve and improve the accuracy reported 
by Fokker-Planck calculations without deep potentials of 
1.7, 0.7 and 0.4 per c ent in energy, mass and angular mo- 
mentum respectively, (|Einsel fc Spurzemlll999. Paper j) . 

There are mainly two effects that make our numerical 
problem much more difficult than earlier ones, increasing 
the errors. First, due to the axisymmetry of the system; 
and secondly, due to the deep growing central potential. It 
made necessary a better resolution in real and velocity space, 
and specially at the time the region of influence of the BH 
dominates the core and during expansion. We increased the 
resolution in the inner parts of the spatial grid (N^,N Z ) 
by setting a linear grid for the core and a logarithmic one 
for the outer regions. The (Nx, AV)-grid has correspond- 
ingly a higher resolution for absolute energies larger than 
E CI it, as mentioned in Section [2TT] We set here Eq = i5 C rit in 
Eq. 1101 A further test of the code was made by reproducing 
the spherically symmetric case with the non-rotating ini- 
tial parameters (as presented in Section I4.1|l . Comparison 
to NBody reali zations are being presented in a forthcoming 
paper. But see I Kim et alJ (|2008l ) for a comparative study of 
non-BH systems. 

As is well known, the long-term evolution of relaxed sys- 
tems does not depend on the details of the initial conditions, 
as these are erased on a relaxation time-scale. To verify this 
we have also started some of our runs from an initial King 



8 J. Fiestas, R. Spurzem 



-- Wo=3.0,w =0.0 

Wo=6.0,w„=0.0 

t -z 




(16.50) 



W„=6.0 



; o.oo) 



Figure 5. Evolution of density at the influence radius for models 
of different initial concentration (Wo = 3.0,6.0). The comparable 
self-similar expansion after collapse shows that the evolution in 
relaxation time scales is not sensitive to the initial configuration 
of the system. Density of model Wo = 6.0 was normalized to its 
density maximum and to its collapse time. Model Wo = 3.0 was 
overplottcd for comparison reasons. 



profile, with a concentration parameter Wo = 3.0 (Model 
M0). The evolution of the density at the influence radius is 
compared to the model Ml (Wo = 6.0) in Fig. [5] Both the 
Wo = 3 and the Wo = 6 models experience a self-similar 
expansion, which approximates n a ~ t~ 2 after collapse is 
prevented. From now on we define the collapse time as the 
time of maximum density at the influence radius in the sys- 
tem (see Tables [3] and [5]) . In Fig. [5] the evolution of model 
Wo = 6.0 was normalized to its density maximum and to 
the collapse time (t cc ), while the model Wo = 3.0 was over- 
plotted for comparison reasons. 

Typical runs for the evolution of one model up to 
~ 50 £ r hi needed about 40 hours on a 3-GHz Pentium IV 
processor (ARI-ZAH, University of Heidelberg). A speed- 
up through parallel processing would be recommended for 
multi-mass versions of the present code (work in progress). 
This performance is not disappointing, taking into account 
that the number of floating-point operations performed per 
time-step in our models is Nx X Ny X x N z and since 
the time steps get much shorter close to t cc . The results 
presented here concentrate on our standard model Wo = 6. 



4 ISOLATED SYSTEMS 
4.1 Spherical symmetry 

In order to test the method, we reproduce the evolution of 
isolated dense stellar systems in the spherically symmetric 
case. We realize this model by setting the initial rotating 
parameter ujo to zero. Mbhj starts growing through accre- 
tion of stars in low-J z orbits. For ~ 0.3£ r M the evolution is 
unaffected by the presence of the small central BH and is 
dominated by the contraction of the core. But the increas- 
ing density supports the growth rate of Mbh in later times, 
when the BH-potential (~ GM^/r) dominates the stellar 
distribution within r a . 

The final steady-state, solid curve in Fig. [6] evolves 



Figure 6. Equatorial density profile (z = 0) for model Ml at 
different times (in parenthesis) given in units of initial half-mass 
relaxation time (t rhi ). Dashed colored lines (black, blue, green 
and red) show the evolutionary profiles and the orange solid line 
the final profile. The dot-dashed line shows the -7/4 slope. The 
location of r a is shown as squares. 



towards a power-law of A = —1.75, according to n a r x . 
This so lution has been exten si vely studied in the spherica l 



1 ins so lution nas been extensivel y studied m tli c spherical 
case bvlBahcall fc WolJ Jl976l) ; lLightman fc Shapiro! (| 19771 ); 
iMarchant fc Shapiro! (|l980D and others. It forms inside r a 
and is maintained in the post-collapse phase, while the evo- 
lution is driven through energy input from the central ob- 
ject, dominating always larger zones. The density profile flat- 
tens close to the center due to the effective loss-cone accre- 
tion and it remains practically unchanged in the halo, where 
the loss-cone loses its significance. 

As the systems evolves, orbits in the region of influ- 
ence of the BH become Keplerian bounded. Their velocity 
dispersion approximates a power-law of -1/2 within the BH 
influence radius r a . Fig. [7] shows the evolution of the total 
one-dimensional velocity dispersion profile in the same way 
as Fig. [6] does for the density. The velocity dispersion grows 
significantly inside r a and faster when the cluster is close to 
collapse, due to the presence of the deep central potential. 
Fig. [8] shows the anisotropy profile in the system at dif- 

2 

ferent times. Anisotropy is defined as A = 2(1 — -^r), where 
the total velocity dispersion of = a 2 > + 2a 2 (o> = a z ). Veloc- 
ity dispersion which is the azimuthal velocity dispersion 
in the direction of rotation and ov are calculated initially 
by taking moments of / with respect to v and v 2 and assur- 
ing conservation of energy and angu lar momen t um p er unit 
volume by encounters between stars (|Goodmanl (|l983l )). The 
initial profile shows a maximum positive halo anisotropy (ra- 
dial orbits dominate the halo in the initial configuration) 
after a very short time (a fraction of trhj- The total amount 
of radial anisotropy, by using the rate K-r/K^, where K r is 
the kinetic energy in the radial degree of freedom and is 
the kinetic energy in the tangential degree of freedom, gives 
a maximum excess of 13 % in K r present in the system at 
the time of the final profile showed in Fig. [8] The small ex- 
cess of K T can be understood, because A{r) rises only in the 
outer regions of the system where the density is low. This 
results are of the same order like the presented in previous 
theoretical studies of anisotropy profiles and evolution of 
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Figure 7. Equatorial profile (z = 0) of ID total velocity disper- 
sion as in Fig.|6]for the density (Model Ml). The dot-dashed line 
shows the -1/2 slope and the dashed lines the evolutionary pro- 
files (as in Fig. [6]l. The location of r a is shown as squares. The 
curve achieves the steady state (orange solid line) in ~ 10t r hi 
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Figure 8. Anisotropy, defined as A = 2(1 — against radius 

for the same model of Figs. [6] and [7] Small tangential anisotropy 
forms inside r a . Radial anisotropy dominates the halo. Curves are 
labeled as in Figs. [6] and [7] 



globu lar clusters |Louis fc Spurzemll99ll ; lGiersz fc Spurzeml 
1 19941 ). Moreover, small negative anisotropy forms slowly in- 
side the BH influence radius (tangential orbits dominate the 
center close to the BH), while rad i al anisotropy remains 
in the halo ilQuinlan G. et all 1 19951 ; iFreitag fc Bend |2002| ; 
iBaumgardt et al.ll20o3 V 



Evolution of Lagrangian radii is a good indicator for the 
contraction and further reexpansion of mass shells. During 
expansion core shells increase as r oc i 2//3 , as expected for 
a system in which the central object has a small mass and 
the energy pr o duction i s con fi ned to a small central volume 
ilHenonl fl965l; IShapird Il977l ; iMcMillan. Lightman fc Cohnl 
Il98ll : iGoodmanl 1 19841 ). In Fig. [9] the influence radius is 
plotted additionally to the Lagrangian radii, rss refers 
to the evaluation of Lagrangian radii at a zenithal an- 
gle, where the effects of probable flattening on the mass 
columns are expected to be less important, that deviations 



Figure 9. Evolution of Lagrangian radii rss, as defined in text, 
for the model Ml. The time is given in units of t r hi. Solid lines 
represent from bottom to top radii of 0.01 %, 0.1 %, 1 %, 5 %, 
10 %, 20 %, 30 %, 50 %, 75 %, 90 % of the total mass of the 
system. The influence radius is additionally plotted (dashed line) 
and the self-consistent expansion is shown by the dot-dashed line 
r oc t 2 / 3 . 



from spherical symmetry are only up to second order in a 
Legrende expansion, i.e. P2(cos9) = 0. That gives 6 = 54, 74 
jEinsel fe Spurzem|[l999. Paper J) . 

Figs (3] and [4] show how Mbh reaches a nearly constant 
fraction of M c n at collapse time (t cc ), while the star accre- 
tion rate (dM/dt) is maximal at t cc due to the higher den- 
sity of orbits in the core decreasing afterwards very rapidly 
(Figs. [3p and[3p). For a density power-law of A = —1.75, 
the expect ed proportionality dM/dt oc t a turns ou t to be 
a = —1.2 (|Amaro-Seoane. Freitag fc Spurzemll2004l ) . 

During evolution, the core is heated via the consump- 
tion of stars in bound, high energetic orbits in the cusp. 
Energy flux is achieved by small-angle, two-body encoun- 
ters, by which some stars lose energy and move closer to the 
BH being eventually consumed, while the stars with which 
they interact gain energy and move outward from the cusp 
into the ambient core. Angular momentum transport is ini- 
tially enhanced by gravo-gyro instabilities and not affected 
by BH accretion of stars on orbits of low J z . Later, when 
core density grows, mass growth rate increases strongly due 
to core contraction to higher densities and stronger stellar 
interaction. 

The general behavior confirms previous studies of spher- 
ically symmetric systems. iMarchant fc Shapiro! (|l98 0l follow 
the evolution of a star cluster containing a central BH (in- 
cluded in their simulations at collapse time). The BH mass 
stalls after approximately 2 relaxation time units to a final 
mass of ~ 4OOOA/0. In our models, a similar rapid evolu- 
tion before expansion is observed, and the final masses are 
comparable (see Tables [5] and [SJ) . 



4.2 Axisymmetric isolated systems 

Evolution of the density profile of model M4 is shown in 
Fig. 1 101 The extent of the BH gravitational influence is 
marked on each curve at the position of r a (squares). Evo- 
lutionary profiles are represented by dashed curves. Note 
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Figure 10. Equatorial density profile (z = 0) for model M4 
(Wo = 6.0, wo = 0.9) for different times given in units of f r hi. The 
orange solid line shows the final profiles and the dashed lines the 
evolutionary profiles. The dot-dashed line shows the -7/4 slope. 
Colors are like in Fig. [6] The location of r a is shown as squares. 



that the limit between the cusp and the core is located at 
~ r a . The central cusp in the density profile grows first very 
slow and faster towards core collapse. It approaches the -7/4 
cusp, like in the spherically symmetric case. 

Fig. [TT] shows the evolution of density contours in the 
meridional plane (zu,z). In the regions where BH star ac- 
cretion dominates (i.e. inside r a ), the isodensity contours 
grow stronger due to the presence of the BH (lighter zones 
in Fig. 111) 1. The cusp forms a strong gradient towards the 
center, forming very fast close to core collapse. Note that the 
flattened shape of the system remains at later times, dur- 
ing the expansion phase (Fig. I11H ) . The velocity dispersion 
is Keplerian within the BH influence radius r a , as Fig. 1121 
shows. Extension of r a is comparable to non-rotating mod- 
els but cusp formation time is shorter, the higher the initial 
rotation parameter (from comparison to Fig[7]) due to the 
faster evolution of this models. Moreover, the rate K T /K t is 
larger the higher the rotation, with maximum values of 1.18, 
1.21, 1.23 and 1.26 for the models M2 to M5, respectively. 

Evolution of Lagrangian radii containing the indicated 
fractions of the initial mass is shown in Fig. [T3]in comparison 
to the non-rotating model. Lagrangian radii give also a qual- 
itative description of the interaction of a growing BH and the 
cluster mass shells. Initially, the BH mass growth is slow due 
to the low central density, and Lagrangian radii are domi- 
nated by core contraction. Finally the collapse is halted and 
reversed and the mass shells re-expand. It happens faster for 
higher rotating models. Note that the smallest radius con- 
tains only 0.01 % of the cluster mass, in order to follow the 
evolution of mass shells closer to the BH. Our single mass 
rotating models show some deviations of the self-similar ex- 
pansion phase, which should be further investigated. 

Fig. 1 141 shows the evolution of density at r a for models 
with rotation parameters ojo=0.0 (non-rotating), 0.3, 0.6, 0.9 
and 1.2. After a similar initial evolution, collapse is faster the 
higher the initial rotation. Gravothermal and gravo-gyro in- 
stabilities drive collapse, while angular momentum is trans- 
ported out of the core in a progressively more efficient way 
for the higher rotating models. Both instabilities occur to- 
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Figure 11. Evolution of density distribution in the meridional 
plane for a model M4. Cylindrical coordinates (to, z) are used. 
Lighter zones represent higher isodensity contours. The time is 
given in units of initial half-mass relaxation time (i r hi). 
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Figure 12. Equatorial profile (z = 0) of ID total velocity dis- 
persion as in Fig. UOI for the density (Model M4). The dot-dashed 
line shows the -1/2 slope and the dashed lines the evolutionary 
profiles. The location of r a is shown as squares. 



gether and support each other. The collapse is reversed due 
to the energy source built by the star accreting BH, while 
the central density drops during expansion. Because radial 
anisotropy first dominates (as shown in Fig. [S| , it supports 
BH accretion of stars in the core, which are able to interact 
with these low-J z (eccentric) orbits in the outer parts. 

As seen in Table [3] collapse times for non-BH models 
are comparable to the BH rotating models and are as well 
shorter for higher initial rotation. t cc varies from 12.20t r h 
(non-rotating models) to 4.6t r h (high rotating models). At 
t C c angular momentum diffusion is more effective, due to 



Dynamical evolution of rotating dense stellar systems with embedded black holes 11 



W„=6.0 



a =0.9 

- «o=1.2 




Figure 13. Evolution of mass shells (Lagrangian radii rss) for 
the models Ml, M4 and M5. Solid lines represent from bottom to 
top radii of 0.01 %, 0.1 %, 1 %, 5 %, 10 %, 20 %, 30 %, 50 %, 
75 %, 90 % of the total mass. Self-similar expansion is shown by 
the dotted line. The time is given in units of t r ^. 
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Figure 14. Evolution of density at r a . Density maxima from right 
to left correspond to the higher rotational parameter uj. The time 
is given in units of t rm . 



the interplay between dynamical instabilities and BH star 
accretion. 

Table|3]shows the final BH mass (in units of Mq) of each 
model and the respective maximal accretion rates (in solar 
mass per year) at the time Mbh stalls and the accretion rate 
begins to slowdown. For a system of M ch = 5- 10 6 M Q , M{f h a11 
varies between 7.5 • 10 3 Mq and 1.5 • 10 4 Mq, which agrees 
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Table 3. Comparison of collapse times t c 
and non-BH models 
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Table 4. Evolution of mass parameters. Column 1: Model name; 
Column 2: Final Mbh in units of initial mass as it reaches an 
asymptotic mass; Column 3: M bh in M0-units; Column 4: max- 
imal accretion rate dM/dt in units of initial mass and half-mass 
relaxation time; Column4: dM/dt in (MQ/yr)-units. 



with the IMBH estimated by theoretical studies and ob- 
servations of globul ar clusters (jGebhardt et alj |2000| . |2002| ; 
iGerssen et al. |2002| ). as expected according to the initial 
conditions of our models (Sect. 13.1(1 . Physical units were 
derivated as described in 13.11 using the initial parameters of 
the corresponding King models (Table [TJ). The general be- 
havior exhibits a slightly decreasing M^ 1 but higher mass 
growth rates for higher initial rotation. In rotating models, 
a stalling of Mbh happens always faster the higher uio is 
(Table El) 

It is known, that in rotating models without BH, the 
total collapse time is shortened by the gravo-gyro effect 
(|Hachisull 19791 . [l982h . by which large amounts of initial rota- 
tion drive the system into a phase of strong mass-loss while 
it contracts (the core rotates faster although angular mo- 
mentum is transported outwards). At the same time, the 
core is heating, while the source of the so-called 'gravo- 
gyro' catastrophe is consumed and the growth in central 
rotation levels off after 2 - 3 t r h towards core collapse 
l|Einsel fc Spurzeml Il999. Paper j) . Simulations into post- 
col lapse phase , driven by thr ee-body binary heating shown 
bv iKim et al] (|2002 Paper ij ) exhibit a faster evolution for 
rotating models. BH rotating models experience in a similar 
way, the onset of gravo-gyro instabilities, as angular mo- 
mentum diffuses outwa r ds, le ading to an increase of central 
rotation l|Hachisulll979l . ll982h . Moreover, BH mass growing 
causes the expansion of the system (Fig. [9] and Fig. 1130 . and 
leads to an ordered motion of high- J z bounded orbits around 
the central BH (tangential anisotropy) supporting develop- 
ment of central rotation. Nonetheless, as the BH reaches its 
final mass, angular momentum continues being transported 
out of the core. Fig. [15] shows snapshots of the evolution 
of the 2-dimensional distribution of v ro t in the meridional 
plane, at representative times, where the lighter areas rep- 
resent contours of higher rotation. Note that an important 
amount of central rotation is still present during the time of 
expansion Fig. I15H (and compare to Fig. I23[l . 

The rate of rotational velocity over velocity dispersion 
represents the importance of ordered motion in comparison 
to random motion. In Fig.[l6]the evolution of the rate V TO t/cr 
at the Lagrangian radii is shown. Up to collapse, there is no 
considerable influence of the BH, while during the expansion 
Kot/o" grows slightly in time, specially for the inner shells. 

The results presented here support the thesis that the 
formation of a massive central dark object, could predict the 
remaining of central rotation in GCs over long evolutionary 
time scales. Nonetheless, the central V IC ,t/& finally decreases 
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Figure 15. 2D rotational velocity distribution in the meridional 
plane for the model M4 (6.0,0.9). 
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Figure 16. V vo tlo at the Lagrangian radii as a function of time 
for the model M4 (6.0,0.9). Vrot/cr at radii of 0.01 %, 0.1 %, 
1 %, 5 %, 10 %, 20 %, 30 %, 50 %, 75 % of the total mass are 
represented by solid lines from bottom to top. The time is given 
in units of t^. 

due to the growing central velocity dispersion, and falls later 
after t cc together with V TOt , while angular momentum is car- 
ried away from the system. In the outer regions the effect is 
smaller, maintaining a slower decreasing. 



5 TIDALLY LIMITED MODELS 

In this models, mass-loss is included, allowing the escape of 
stars through the energy tidal limit (see Fig. [T} . While Mbh 
grows and central density increases within r a , the system 
loses mass through the outer tidal boundary due to relax- 
ation effects. Evolution of the density profile of model M4T 
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Figure 17. Equatorial density profile (z = 0) for model M4T 
(tidally limited case). Times are given in parenthesis and in units 
of t rhi . The dot-dashed line shows the -7/4 slope and the dashed 
lines the evolutionary profiles. Solid (orange) line shows the final 
profile. The location of r a is shown as squares. 



is shown in Fig. [17] The extent of the BH gravitational influ- 
ence is marked on each curve at the position of the influence 
radii r a (squares). Evolutionary profiles are represented by 
dashed curves. The power-law of A = —1.75 is not com- 
pletely reached as in the isolated case (Fig. I17p . probably 
due to the strong mass loss. In the halo, it gets steeper be- 
yond r a up to the tidal radius (rud), where the loss-cone 
loses its significance, rua itself becomes smaller in time as a 
consequence of tidal mass-loss. 

Fig. [18] shows the evolution of density in the merid- 
ional plane (zu,z). In the regions where BH star accretion 
dominates (i.e. inside r a ), the isodensity contours grow due 
to the presence of the BH (lighter zones in Fig. I18[) . Note 
that scales are different for the bottom figures due to the 
shrinking of the system which loses mass through the outer 
tidal boundary, and the cluster tidal radius becomes smaller 
(darker areas) . The shape of the system becomes faster more 
spherical that in the isolated case. Like in the isolated case, 
orbits in the region of influence of the BH become Keplerian 
bounded. Their velocity dispersion grows towards a power- 
law of -1/2 within the BH influence radius r a (Fig. 1 19)1 . 

Fig.[20]shows the anisotropy profile in the system at dif- 
ferent times. The initial profile shows radial halo anisotropy 
(radial orbits dominate the halo in the initial configuration). 
Tangential anisotropy seems not to form inside the BH in- 
fluence radius, as it does in the isolated case. Moreover, at 
later times, tangential orbits dominate the halo as a con- 
sequence of an effective J z -transport outwards and the ac- 
cretion of preferentially radial orbits by the central BH. In 
general, a faster evolution in higher rotating models leads 
to a smaller M^ff 11 (see Table [6]), and thus to a smaller tan- 
gential anisotropy (i.e. smaller in model M4T than MIT). 
This was expected, since a more massive BH will consume 
more stars in preferentially radial orbits and the life time of 
high rotating tidally limited systems is too short to develop 
higher Mbh- Note that at this evolutionary times (t ~ 5t r hi 
for M4T) the cluster has lost more than 50 % of its mass 
(see Table [5} and the densities in the outer parts are much 
lower than in the isolated models. As a consequence, the 
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Figure 18. Evolution of density distribution in the meridional 
plane for a model M4T in the tidally limited case. Cylindrical 
coordinates (zu, z) are used. Lighter zones represent higher iso- 
density contours. Note that scales are different in the bottom 
figures due to the shrinking of the outer tidal radius. 
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Figure 19. Equatorial profile (z = 0) of ID total velocity dis- 
persion (Model M4T) in the tidally limited case. The dot-dashed 
line shows the -1/2 slope. Evolutionary profiles are labeled like in 
Fig. 1171 The location of r a is shown as squares. 



measured rates Kr/K^ range from 0.98 to almost 1.00 for 
all these models. At the same time, r a becomes larger and 
dominates almost the hole system, which itself is close to 
dissolution (see Fig. [2Tj) 

Fig. [21] shows the evolution of central density for all 
models, with rotation parameters cjo— 0.0 (non-rotating), 
0.3, 0.6, 0.9 and 1.2 for the tidal limited case. Collapse time 
is reached faster, the higher the initial rotation, in compari- 
son to the isolated models. Gravothermal and gravo-gyro in- 
stabilities support each other and collapse phase is reversed 
due to the energy source built by the star accreting BH. 



W„=6.0 



Figure 20. Anisotropy, defined as A = 2(1 — -^r), against radius 
for the same model of Figs. [T71 and I19l Dashed curves are evolu- 
tionary profiles (labeled with times in parenthesis). Solid (orange) 
curve (latest evolutionary time) is labeled separately. 
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Figure 21. Evolution of density at the influence radius for models 
of Wo = 6.0 and initial rotation ujq =0.0,0.3,0.6,0.9,1.2 (density 
maxima from right to left). 



Evolution of Lagrangian radii containing the indicated 
fractions of the initial mass is shown in Fig. 1221 Due to mass- 
loss, the outer mass shells are rapidly truncated, the faster 
the higher initial rotation. At the same time, the core density 
grows (higher disruption rates). Later, collapse is halted and 
reversed (while accretion rate slows-down rapidly) and the 
mass shells re-expand for a few t r j,. 

Collapse times for tidal limited models with same initial 
conditions, can be seen in Table[5] where collapse parameters 
of rotating BH and non-BH models are compared. Collapse 
times are comparable to the non-BH models and slightly 
smaller for increasing rotation. The cluster loses, at collapse 
time, between 40 and 85% of its mass in the non-BH models 
and between 48 an 87 % of the initial cluster mass in the 
BH models. In both cases the mass-loss is higher for higher 
rotation. The time at which the cluster loses half of its mass 
is shorter, the higher the initial rotation of the model. The 
effect is driven by the interplay between relaxation effects 
(gravo-gyro instabilities) , BH star accretion and tidal mass- 
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Figure 22. Evolution of mass shells (Lagrangian radii r^) for the 
models MIT, M4T and M5T. Solid lines represent from bottom 
to top radii of 0.01 %, 0.1 %, 1 %, 5 %, 10 %, 20 %, 30 %, 
50 %, 75 %, 90 % of the total mass. Note that some of the outer 
Lagrangian radii in the rotating models disappear at very early 
times, due to strong mass loss. 
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Table 5. Comparison of collapse parameters between rotating 
BH and non-BH models. t C c: core collapse time 
tso : time at which the cluster lost half of its mass 
M cc : current cluster mass at t rj t C c 



loss. This will also involve a faster dissolution of the cluster 
in the galactic tidal field. 

Table [6]shows the final BH mass for each model and the 
respective maximal accretion rates. Ai^f" varies between 
2.3- 1O 3 M and 7.7- 1O 3 M , which are smaller masses than 
in the isolated case. The general behavior exhibits a decreas- 
ing M^f but a higher mass growth rates corresponding to 
higher initial rotation. A stalling of Mt>h is always faster the 
higher uio is. 

The cluster mass-loss due to tidal effects of the par- 





Figure 23. 2D rotational velocity distribution in the meridional 
plane for a model M4T (6.0, 0.9) in the tidally limited case. 



ent galaxy is very strong during the re-expansion of the 
cor e. The acce l eration of mass -loss is similar as the observed 
bv iKim et all (12002 Paper Ilh in the post-collapse models 
driven by binary heating, although the effect in the present 
BH-models is more pronounced, with the consequence of a 
faster evolution of the cluster towards sphericity and final 
dissolution. 

Tidally limited models experience the onset of gravo- 
gyro instabilities, as angular momentum diffuses outwards, 
leading to a strong mass-loss and a experience a limited 
increase of central rotation. At the galaxy tidal boundary, 
mainly circular tidal orbits in the halo are lost, while high 
eccentric (low J z ) orbits can interact with stars in the core. 

Fig. [23] shows snapshots of the evolution of 2- 
dimensional distribution of rotational velocity in the merid- 
ional plane, at representative times, where the lighter ar- 
eas represent contours of higher rotation. Rotation is lost 
stronger than in the isolated model. Angular momentum 
transport and the growing BH-potential support the devel- 
opment of ordered motion in the core and, at the same time, 
trigger mass-loss through the tidal boundary. 
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6 CONCLUSIONS AND OUTLOOK 




Table 6. Description of mass parameters. Column 1: Model 
name; Column 2: Final M^h in units of initial total mass as it 
reaches an asymptotic mass; Column 3: Mbh in M -units; Col- 
umn 4: maximal accretion rate dM/dt in units of initial mass and 
half-mass relaxation time; Column 5: dM/dt in (M /yr)-units. 



The variety of environments in which dense stellar systems 
(GCs, GN) form and evolve, make them target of study of 
fundamental dynamical processes. Observational studies of 
GCs suggest the existence of IMBHs in their centers (like 
Gl and M15) and is well known, that some of them show 
flattening due to system rotation. The improvements in ob- 
servational and theoretical studies of dense stellar systems 
in the last years, has led to a better understanding of central 
BHs and their environments, but at the same time, opened 
new questions due to their complexity. This is the reason, 
why theoretical models are of importance to elucidate the 
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origin of the observed phenomena, to be able to explain their 
formation, evolution and interaction, and predict possible 
evolutionary scenarios, which can be confirmed by observa- 
tions. The presented theoretically formulated evolutionary 
models extend the model complexity of spherically symmet- 
ric systems through the implementation of differential rota- 
tion and BH star accretion. 

Results can be summarized as follows: 

• Core collapse is an evolutionary property in self- 
gravitating s ystems without BH ( Einsel fc Spurzeml 
11999. Paper 7 ! iKim et all 12002 Paper HI) . Gravo-gyro ef- 
fects are coupled to gravothermal instability and drive core 
contraction. We start with a seed central BH, which grows 
over relaxation time scales due to stellar accretion. Evolu- 
tion towards core contraction increases the central density 
and supports star accretion by the central BH. The BH acts 
as an energy source, through which energetic stellar orbits 
are formed, which easily go into the loss-cone or interact 
with stars inside r a being able to reverse collapse, triggering 
core expansion. A rapid ex pansion has been obse rved as 
well in TV-Body realizations (|Baumgardt et aT1l2004l ). which 
use high concentrated initial models (Wo = 10) and higher 
initial BH masses (~ 1 — 5%M c i) in their simulations. They 
show acc retion rates which a gree to the classical approx- 
imation (|Frank fc Reesl |l976j) applied to a Bahcall-Wolf 
cusp, as we also show in Figs. [3] and U In the presence 
of an external potential (tidal limit), a faster evolution 
accelerates mass-loss and leads to cluster dissolution. 

• Final steady state solutions are found in all isolated 
models, which approach the -1.75 slope in the density cusp 
and the -0.5 slope in the velocity dispersion cusp inside the 
BH influence radius r a , corresponding to Keplerian bounded 
orbits, independent of initial rotation. 

• The final Mbh nearly stalls at ~ 0.001M c ii, and grows 
slower during the post-collapse phase, while BH mass ac- 
cretion rate (dM^/dt) decreases strongly after reaching a 
maximum before core expansion, due to the higher density 
of orbits in the core. In the tidally limited model mass-loss is 
very strong during the re-expansion of the core. The cluster 
mass reaches in a fraction of i r hj after collapse values at least 
one order of magnitude smaller than its hosted BH. For a 
cluster of 5 ■ 1O 6 M , M{f h a11 varies between 7.5 • 10 3 M Q and 
1.5 • 10 4 Mq for isolated models of different initial rotation, 
and between 2.3 • 1O 3 M and 7.6 ■ 1O 3 M for tidally limited 
models. As mentioned in Sect. 13. ll this values reproduce the 
expected mass of IMBHs, since we set our initial parameters, 
according to physical properties of these objects. 

• High rotating, moderate concentrated models (M4, M5) 
maintain central rotation at collapse and during the expan- 
sion phase in an efficient way, in comparison with models 
without BH. They are able to maintain an efficient angular 
momentum diffusion, and at the same time are concentrate 
enough to avoid an excessive mass-loss. Both effects support 
the accretion of stars in low-J z orbits. These models show 
a stable evolution of ordered vs. random motion (Vrot/c) in 
BH-models, up to the expansion phase. 

Since flattening supported by rotation is a well known 
phenomena in GCs and observational evidences of the exis- 
tence of central dark objects in some GCs, with and with- 
out rotation, there is a motivation for the study of this 
constraint in the long term evolution. Although some con- 



straints are still missing, like a mass spectrum or stellar 
evolution (work in progress), as well as a more realistic 
criterion for tidal mass-loss ( as observations suggest, e.g. 
iMackev fc van den Berd 12003 ). our models are consistent 
with existent theoretical studies on the general evolution of 
systems embedding BHs and additionally consider the im- 
portance of initial differential rotation, which needs to be 
taken into account for the understanding of GC formation 
and evolution, specially w hen it can be high eno ugh in young 
clusters (e.g. in the LMC. iBrocato et all (|2004h ). 

Moreover, galaxy cores are known to harbor SMBHs 
and some of them are 'collisionaP, in the sense that their 
relaxation times are < 10 10 yr (like M32 and the MW). 
They show cuspy density profiles, which shorten the times 
of relaxation the smaller the distance to the center of the 
nuclei and might support the formation of a steady state 
configuration (Bahcall-Wolf solution) at times of the order 
of 10 Gyrs. 

Since the models presented here have only one mass 
component, the effect of mass segregation in a realistic multi- 
mass system will shorter times of evolution |Gurkan et al.1 
l2004h . leading to a faster set-off of expansion in the presence 
of a central BH. The central velocity dispersion arises be- 
cause of the BH-induced Bahcall-Wolf cusp. This increase 
is not affected by the presence of a spectrum of stellar 
masses but the high mass stars are expected to s how a lower 
cusp in their c entral dispersion (as reported by IKim et all 
12004 Paper lit) , leading possibly to a higher or at least more 
stable Kot/c for this mass classes. Multi-mass models with 
BH are being currently developed and comparison to N- 
Body models are aimed to complement this calculations, 
using the highest particle num ber permitted at the time 
(N ~ 10 6 ) |Berczik et alj|2006h . 

As shown, the amount of rotation present in the system 
during its dynamical evolution is strong influenced by the 
interplay between angular momentum diffusion (gravo-gyro 
instability) and the redistribution of high energy orbits close 
to the BH (loss-cone refilling). Since a central BH is able to 
'consume' angular momentum from the system, in form of 
stars, it might become itself an angular momentum source, 
which could be able to rotate (Kerr Black Hole), permitting 
also a more efficient angular momentum transport outwards, 
through interaction with core stars, driven by relaxation. A 
binary black hole (BBH), could in a similar way, lead to a 
more efficient support in the development of rotation in its 
zone of infl uence, modifying substantially the final shape of 
the cluster (jMapelli et al.ll2005l ; iBerczik et alj|2006ft . 

The models presented, are able to reproduce 2D dis- 
tributions (in the meridional plane) of density, cluster- 
and BH-potential, velocity dispersions, rotational velocity, 
anisotropy, dynamical ellipticity, among other parameters, 
at any time of evolution and deep in the stellar cusp sur- 
rounding the central BH. They make possible the study of 
kinematical and structural parameters in time, which can 
complement and test observational measurements contribut- 
ing to the understanding of the common evolution of star 
clusters and galaxies. 
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